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Abstract 

We report a first-principles description of the induced wake potential and 
density in real solids. The linear-response formalism is used to obtain the 
potential and density induced by an external charge penetrating through an 
inhomogeneous periodic system. The linear dynamical response of the system 
is evaluated in the random-phase approximation, by including the full band 
structure of the solid. The impact of the periodic crystalline potential is 
analyzed, and variations of the wake along different channels in Al and Si are 
investigated. 
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I. INTRODUCTION 



A charged particle penetrating a solid causes a distortion in the electronic density around 
the particle and behind its position: this is what Bohr called the induced density wake [|l|. 
Related to this induced density there is an induced potential. For a sufficiently high velocity 
the induced wake shows an oscillatory behaviour. When the velocity of the projectile is 
larger than the average velocity of target electrons (typically vp-, the Fermi velocity), one 
may consider a linear response of the medium. However, in the case of projectiles moving 
with smaller velocities, nonlinearities may play an important role for metallic densities (2 < 
Ts < 6, Ts being the average electron distance 0). 

Much research has been oriented to the study of these quantities. Pioneering work on 
dynamic screening was performed by Neufeld and Ritchie 0. They evaluated the induced 
potential and density wake by using a local dielectric function as the linear response function 
of the medium. Based upon a linear response of the target, different aspects related to the 
induced polarization such as wake riding states or the spatial distribution of the induced 
potential and density have been investigated. The induced potential and density have 
been studied beyond linear-response theory, in the static electron gas approximation |§] 
and within hydrodynamical formulations [p|-p!T[]. The quadratic induced polarization by an 
external charge in the full random-phase approximation (RPA) has been recently reported 
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In previous works the induced potential and density have been evaluated on the basis 
of a jellium model of the target, which consists of an isotropic homogeneous electron gas 
embedded in a uniformly distributed positive background. However, in a more realistic 
approach valence electrons move in a periodic potential and one-electron excitations split 
into the so-called energy bands. The impact of band-structure effects on plasmon dispersion 
p!4|-[T^, dynamical structure factors ||T7|-pO[|, stopping power [|21[, and hot-electron 



curves 



lifetimes p2| , |23[| has been investigated only very recently. It has been shown that these 
effects can be important even in the case of free-electron metals such as aluminum. The 
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wake potential induced by swift protons through different sohds has been studied using a 
model dielectric function ||24|| . 

In this paper we report, within linear-response theory, an ab initio evaluation of the 
induced potential and density of ions moving through real solids. Within our description 
we can obtain the induced potential under channeling conditions, i. e., when the projectile 
penetrates through different symmetry directions of the solid. In section II we derive ex- 
plicit expressions for the position-dependent induced potential and density. In section III 
numerical calculations of the induced potential in Al and Si are presented, for both random 
and crystal-symmetry incident directions. In section IV the most relevant conclusions are 
summarized. Atomic units are used throughout, i. e., me = e = ^ = 1. 

II. THEORY 

We consider a point particle of charge Zi moving in an inhomogeneous system with 
velocity v and impact vector b, such that 

p""*(r,t) = Zi5(r-b- vt). (1) 

For a periodic crystal we can write: 

pg*(q, uj) = 27rZie'^('i+^)-^5 [cj - (q + G) ■ v] , (2) 

where G is a reciprocal-lattice vector and q lies in the Brillouin zone (BZ). The external 
charge induces a density, so that the total density variation of the medium is given by the 
sum 

PG(q, uj) = p^cf*(q, Lo) + p'S^(q, u). (3) 
Poisson's equation allows us to write the potential related to a given density. In our case, 

0G(q,w) = t;G(q)PG(q,^), (4) 



3 



where fG(q) = 47r/|q + Gp and 0G(q, i^) represent the Fourier components of the bare 
Coulomb potential and the total potential, respectively. Within linear response theory, we 
can write the total potential in terms of the external charge: 

0G(q,^) = 5IeG!G'(q>^)^G'(q)PG*(q>^)- (5) 

G' 

Cq^q/ (q, to) are the Fourier coefficients of the dielectric function, which within the RPA read 

eG,G'(q,^) = f^G.G' -t'G'(q)XG,G'(q>^)> (6) 

where Xg c (q> ^) represent the Fourier coefficients of the density- response function of non- 
interacting electrons, 

/ \ 1 \ ^ \ ^ \ ^ /k,n /k+q,n' 

XG,G'(q,^^) =nl^2^2^ 



^ BZ n n' ^k,n " £^k+q,n' + (^^ + 
x(0k,n|e-^(''+^>'-|0k+q,n')(0k+q,n'|e'(''+^')"-|0k,n). (7) 

Here, the sums run over the band structure for each wave vector k in the first BZ, |0k,n) 
and ek,n are the one-electron wave functions and energies, and /k,r), are Fermi factors, /k,n = 
6{Ef — ^k.n)- ^ is the normalization volume. 

Combining Eqs. (0), and we obtain the following equation for p*""^ in terms of 

p-'^(q,.;) = [|q + G|2eG;c;,(q,c.)|q+ G'r' - ^G.G'] PG'(q,^)- (») 
G' 

Substituting Eq. (^) into Eq. (|^) and Fourier transforming back to real space yields 

7 BZ 

^ q G K 

|q + G|'eG!G+K(q, (q + G + K) ■ v) |q + G + - 5g,g+k] , (9) 
where we have made use of 



[3q_^^^^ 
(27r)3 ~ Q 



''"-^^Es:. (10) 
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and we have set G' = G + K. 

If we consider a definite trajectory of the projectile, only the K vectors such that K- v = 
contribute to the sum, and 
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(11) 



|q+G + K|2 |q + G|2_ 

The prime in the sum over K accounts for the K ■ v = condition. Therefore, Eq. ([llD gives 
the density induced by an external charged particle moving in a definite trajectory through 
a crystal. 

In order to get the induced potential we just apply Poisson equation, which yields 
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eG!G+K(q> (q + G) ■ v) 
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(12) 



|q+G + K|2 |q+G|2_ 

We can define an average or induced random potential as the mean value over impact vectors 
b, which results in the K = term of Eq. (P^): 

-4nZi gi(q+G)-(r-vt) 



EE 



eG!G(q, (q+G)-v)-l 



(13) 



q G lq + Gp 

The most important contribution to the position-dependent induced potential of Eq. {\L2^ 
is provided by the K = term. For those directions for which there are not reciprocal vectors 
satisfying the K ■ v = condition, we have the average potential of Eq. (|1^). For a few 
highly symmetric or channeling conditions, non-negligible corrections to the random result 
are found. The random induced potential exactly coincides with the well-known jellium 
result 1^ when one replaces the inverse dielectric matrix entering Eq. (|T3D by the inverse 



dielectric function of a homogeneous electron gas ||2^ . 

The position-dependent stopping power can be obtained directly from Eq. (|T^), by 
simply taking into account that 



-f-fv^'"t-«-v, (14) 

which gives an expression that exactly coincides with the position-dependent stopping power 
derived from the knowledge of the imaginary part of the projectile self-energy ||2l|| . 

The main ingredient of our calculation of the induced potential and density is the inverse 
dielectric matrix, which we have evaluated in the RPA by inverting Eq. @. Hence, at this 
point we have only considered the average electrostatic interaction between the electrons. 
We have not found differences in the induced potential by including many-body short- 
range correlations in the response of the system. The one-electron Bloch states entering 
Eq. d^) are the self-consistent LDA eigenf unctions of the Kohn-Sham equation of density- 



functional theory (DFT) p7| , |28| . We first expand the states in a plane- wave basis and 
then solve for the coefficients of the expansion self-consistently. The electron-ion interaction 
is described in terms of a non-local, norm-conserving ionic pseudopotential [^]. The XC 
potential is computed with use of the energy first calculated by Ceperley and Adler and 
then parametrized by Perdew and Zunger We subsequently evaluate the Xg ci^^^) 
polarizability and invert Eq. as we sum over q, G, and K to obtain the induced potential 
and density. 



III. RESULTS 



In this section we present the results of our calculations of the induced potential when a 
charged particle penetrates through aluminum and silicon. The average electron densities of 
Al and Si are similar: the corresponding free-electron gas (PEG) is characterized by = 2.07 
for Al and = 2.01 for Si. However, aluminum is a metal and silicon a semiconductor and 
the induced potential will exhibit a different behaviour, which we would not obtain on the 
basis of a PEG calculation. 

Although Al is usually regarded as a jellium-like material, inelastic X-ray scattering 
experiments P2| and theoretical analyses of the dynamical structure factor |p!7| , p!^ , stopping 
power [^, and hot-electron decay rates |^ have revealed that it is necessary to take into 
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account the full band structure for a proper description of its electronic properties. In this 
work, Bloch states have been expanded in a plane-wave basis with an energy cutoff of 12 Ry, 
which corresponds to keeping in this expansion ~ 100 plane waves. The sums over the BZ 
for both the polarizability Xgg'(^''^) ^^"^ induced potential have been performed on 
10 X 10 X 10 Monkhrost and Pack meshes [^. The sum over reciprocal-lattice vectors in the 



potential has been extended to the first 15 G vectors, which corresponds to a cutoff in the 
momentum transfer of 2.9qF {qf is the Fermi momentum). We have included 30 bands in 
the sums over the band structure for each k vector in Eq. (^, which allows us to calculate 



the induced potential and the stopping power up to velocities of the order of 2 a.u. . 

Silicon is a covalent crystal which shows strong valence-electron density variations in 
certain directions. This allows the formation of channels in which the density is very low. 
For example, the integrated density in the (110) channel varies from = 1.49 at the atomic 
row to Ts = 3.37 at the center of the channel, i.e., an 80%. This has an obvious impact on 
the induced potential and the stopping power. The covalent character of Si imposes a higher 
cutoff for the Bloch-state expansion than in Al, due to the higher degree of localization of 
the electronic states in the former material. This also results in Si having bands that are 
flatter than in the case of Al, so that the number of bands included in the calculation of 
the polarizability Xgg'(^''^) niust be larger in the case of Si if the same energy transfers 
as in Al are to be included. We have used a cutoff of 16 Ry (~ 300 plane waves/state), and 
have included 100 bands in the calculation of Xg g'(^' "^^^ sums over the BZ have been 
performed on 8 x 8 x 8 Monkhorst and Pack meshes ||3^ , and the sum over reciprocal-lattice 
vectors has been extended to the first 15 G vectors, which corresponds to a momentum- 
transfer cutoff of 2.1qF- 

We focus on the spatial distribution of the induced potential along the incoming particle 
trajectory. The z coordinate appearing in the figures is always relative to the particle 
position. The stopping power derived from the slope of the induced potential at the projectile 
position coincides with the stopping power reported in Refs. |]2T]] and . An ion of Zi = 1 
is always considered. 



In Fig. 1 we have plotted the random potential induced by an ion with v = 0.6 a.u. 
through Al (a) and Si (b). The solid line represents the calculation for the real solid, 
whereas the dashed line represents the FEG calculation. For this low-velocity regime (below 
vp), differences between full band-structure and FEG calculations come from the presence 
of interband transitions and, also, the gap in the case of Si. The splitting of the bands makes 
the polarization easier and, for this reason, the induced potential at the origin is higher for 
the real crystal than for the FEG in Al. In the case of Si, the slope of the induced potential 
at the origin is smaller than in the corresponding FEG, due to the presence of the gap, 
which, for low-energy transfers diminishes the polarization of the electron system. The gap 



of semiconductors like Si brings about a non-linear stopping power for low velocities |34 
and an induced potential which presents, at these velocities, a lower slope. 

In Fig. 2 we have plotted the random potential induced by an ion with v = 1.5 a.u. 
through Al (a) and with v = 1.6 a.u. through Si (b). Solid and dashed lines represent 
full band-structure and FEG calculations, respectively. These velocities are well over the 
plasmon-excitation threshold {vt ~ 1.3 a.u. in both cases). An oscillatory behaviour appears 
behind the ion due to plasmon excitation [§[^0], the wavelength of the oscillations being 
of ~ 2TTv/ujp {ujp is the plasma frequency). As can be appreciated in Figs. 2a and 2b, the 
wavelength of the random potential is the same in the real crystal and in the corresponding 
FEG for both materials. This is a consequence of the fact that plasmon contributions to 
the stopping power are the same for both the real solid and the FEG [^]. However, the 
oscillations in real Si are more damped than in real Al, due to the shorter lifetime of plasmons 
in Si which stems from the higher density of bands, thus increasing the decay channels . 



In Fig. 3 we have plotted the position-dependent potential, Eq. (jT^), induced by an 
ion with v = 1.5 a.u. along the (100) direction in Al (a), and with v = 2.0 a.u. along 
the (110) direction in Si (b). The impact vector is 6(0,1,0) for Al and 6(— 1,1,0) for Si. 
b is measured in units of the lattice constant ac of the target. Calculations for an impact 
parameter 6 = (atomic row) and b = 0.25 (center of the channel) are represented by solid 
and short-dashed lines, respectively. For comparison, a local-density approximation (LDA) 
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35| of the induced potential is also displayed in Fig. 3 by a dashed line for 6 = and 



a dotted line for b = 0.25. In this approach the position-dependent induced potential is 
obtained as the potential induced in a FEG with an electron density equal to the average 
electron density along the projectile path. The slope of the induced potential gives the 
stopping power in each case. We can appreciate a substantial variation in the slope of Si as 
we move across the channel, due to the great variation of the density. A lower density at 
the center of the channel results in a smaller slope and a lower stopping power than in the 
atomic row and in the random case. Variations in the spatial distribution of the induced 
potential in Al are fairly small as b changes, because the electron density is almost flat 
in this material. If we calculated the potential at the ion position as a function of b, we 
would obtain a qualitatively similar result in both the ab initio and FEG evaluation, i.e., 
following the local density variation |21|. However, it is clear from Figs 3a and 3b that the 
spatial distribution of the induced potential is not well described by the LDA calculation. 
Above all in Si, Fig. 3b, the differences are more dramatic. The oscillations behind the 
ion have the same wavelength independently of b in the crystal calculation, whereas in the 
LDA calculation the wavelength depends on the average density for each path. Plasmon 
contributions appear in the K = term of the position-dependent induced potential of 
Eq. ([T^) . The remaining terms (K ^ 0) are mainly corrections to the single particle e-h 
contribution, which are due to the presence of density variations away from the projectile 
path, i.e., the so-called crystalline local- field corrections. For this reason, the wavelength of 
the induced potential remains unchanged. Furthermore, the stopping-power peak is located 
at the same position for any velocity direction and for any impact parameter [p^ 



smce 



plasmons are collective excitations which involve all the electrons of the system. However, 
in an LDA calculation the stopping-power peak is located at the velocity that corresponds 
to the average density of the path and not to the average target density. These results lead 
us to the conclusion that LDA calculations are not suitable for the calculation, at these 
velocities, of the position-dependent induced potential and stopping power. 
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IV. CONCLUSIONS 



We have presented full band-structure calculations of both random and position- 
dependent induced potentials in Al and Si. The linear-response formalism has been used 
to obtain the potential and density induced by an external charge penetrating through an 
inhomogeneous periodic system. The random potential has been evaluated in the RPA for 
velocities below and above the plasmon-excitation threshold. For low velocities, differences 
between the FEG and ab initio calculations come from the sensitiveness to the band struc- 
ture of the target. For higher velocities, we have shown that oscillations behind the ion have 
the same wavelength in both the FEG and the real crystal, due to the fact that plasmon ex- 
citation remains unchanged. Finally, we have investigated the position-dependent potential 
induced by projectiles incident along the (100) direction in Al and along the widest channel 
in Si, the (110) direction. Variations in the spatial distribution of the induced potential 
are more pronounced in the case of Si. Besides, we have shown that the LDA calculation 
does not properly account, at the velocities under study, for the spatial distribution of the 
induced potential along the channel. 
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FIGURES 

FIG. 1. Random potential induced by a proton {Zi = 1) moving with v = 0.6 a.u. through 
Al (a), and Si (b). Solid and dashed lines represent full band-structure and FEG calculations, 
respectively. 

FIG. 2. Random potential induced by a proton (Zi = 1) moving with v = 1.5 a.u. through 
Al (a), and with v = 1.6 a.u. through Si (b). Solid and dashed lines represent full band-structure 
and FEG calculations, respectively. 

FIG. 3. Potential induced by a proton [Zi = 1) moving with v = 1.5 a.u. along the (100) 
direction in Al at the impact vector b = 6(0,1,0) (a), and with v = 2.0 a.u. along the (110) 
direction in Si at the impact vector b = 6(— 1, 1, 0) (b). Solid and short-dashed lines represent full 
band-structure calculations for 6 = and b = 0.25, respectively. Dashed and dotted lines represent 
LDA calculations for 6 = and b = 0.25, respectively. 
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Figure 1 (a) 




Figure 2 (b) 
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